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GENERAL  DISCUSSION 


Science  Applications,  Inc.  has  utilized  the  family  of 
Flux-Corrected  Transport  (FCT)  algorithms  at  NRL  to  model 
the  airblast  produced  by  nuclear  and  non-nuclear  events. 
The  FCT  algorithms  combine  low  order  numerical  integration 
schemes  which  are  non-dispersive  with  high  order  schemes 
that  are  non  diffusive.  The  net  effect  is  to 
significantly  improve  the  solution  of  flow  fields  with 
sharp  gradients. 

The  one-dimensional  code  can  be  configurated  for 
cylindrical,  spherical  or  planar  geometries.  This  code  is 
primarily  used  to  generate  the  flow  field  from  spherical 
high  explosive  charges.  The  two-dimensional  code  is  a 
cartesian  or  cylindrical  code.  Both  codes  use  the  same 
FCT  algorithm.  The  three-dimensional  code,  however, 
employs  a  leap  frog  FCT  scheme.  It  is  not  as  accurate  but 
runs  a  factor  of  two  faster  on  a  computer. 

One  and  two-dimensional  simulations  of  PBX  9404  were 
completed  and  comparisons  made  with  data.  Appendix  A 
details  the  results.  A  two  dimensional  calculation 
approximating  an  infinite  plane  of  explosions  was 
completed  and  is  reported  in  Appendix  B.  The  limitations 
and  comparisons  of  both  nuclear  and  non-nuclear  events  are 
outlined  in  Appendix  C.  When  blast  waves  occur  over  dusty 
surfaces  there  is  significant  transport  of  the  dust. 


Appendix  D  reports  the  results  of  research  into  the  dust 
and  shock  interaction.  Finally,  an  attempt  was  made  to 
approximate  the  gas  dynamic  equations  through  power  series 
expansions  for  the  case  where  a  planar  shock  encounters  a 
wedge.  Appendix  E  explains  the  analysis. 
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ADAPTATION  OF  FLUX-CORRECTED  TRANSPORT  ALGORITHMS 
FOR  MODELLING  BLAST  WAVES 


Blast  wave  phenomena  include  reactive  and  tvo-phase  flows  associated 
with  the  motion  of  chemical  explosion  products;  the  propagation  of  shocks, 
rarefaction  waves,  and  contact  discontinuities  through  a  nonideal  medium 
(real  air,  possibly  thermally  stratified  and  containing  dust  and  water 
vapor);  and  the  interaction  of  the  blast  waves  (including  boundary  layer 
effects)  with  structural  surfaces.  Flux-Corrected  Transport  (FCT)  represents 
ar.  accurate  and  flexible  class  of  methods  for  solving  such  nonsteady 
compressible  flow  problems  (Boris  and  Book,  1976).  Coupled  with  a 
nondiffusive  adaptive  gridding  scheme  (Book,  et  al. ,  19S0;  Fry,  et  al. , 

1981) ,  it  enables  complex  time-dependent  shocks  to  be  efficiently 
"captured.” 

In  models  which  treat  all  the  physical  effects  required  for  blast  wave 
simulation,  truncation  errors  inherent  in  the  underlying  finitei-difference 
scheme  are  exacerbated  by  nonlinear  coupling  between  the  fluid  equations  and 
by  the  greater  complexity  of  the  phenomena  being  simulated.  Typical  of  these 
errors  are  the  "terraces"  which  develop  under  some  circumstances  on  the 
flanks  of  sloping  profiles  when  the  growth  of  ripples  due  to  phase  errors  at 
short  length  scales  is  terminated  by  the  action  of  the  flux  limiter.  Two 
approaches  are  possible  toward  eliminating  them:  improving  the  short- 
wavelength  phase  and  amplitude  properties  of  the  underlying  algorithm,  and 
switching  on  additional  diffusion  locally.  The  latter  approach  folds 
information  about  the  shape  of  the  profile  and  the  nature  of  the  physical 
process  taking  place  (e.g. ,  rarefaction)  into  the  switch  criterion,  thus 
changing  the  FCT  technique  from  a  "convective  equation  solver"  to  a  "fluid 
system  solver."  In  doing  this,  care  mist  be  taken  to  avoid  losing  the 
accuracy,  robustness  and  problem-independence  which  constitute  valuable 
attributes  of  FCT  algorithms  (Book,  et  al. ,  1981). 

Tests  carried  out  on  scalar  advection  of  simple  density  profiles  by  a 
uniform  flow  field  show  that  terracing  does  not  require  either  diverging 
velocities  or  discontinuities  in  the  profile,  but  appears  typically  (for  v  > 
0)  where  the  first  and  second  derivatives  of  density  have  the  same  sign 
(Fig.l).  In  order  to  improve  the  properties  of  the  basic  difference  scheme, 
we  propose  a  new  algorithm  for  integrating  generalized  continuity  equations 
over  a  timestep  6t.  Consider  the  following  three-point  transport  scheme: 

V  pr  n(p3+i  •  p°-i]  +  k(p°+i  ■  2p3  -  p3-i); 

Pj  *  3^  -  e(p°+1  -  P°_i)  +  *(p°+i  -  +  pj-i); 

p3  =  -  y(*J+l/2  ‘  *J-l/2 

where 

Vl/2  =  PJ+1  “  V 
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The  arrays  {p®}  and  {p1?}  are  the  old  and  nev  densities,  Pj  and  p<  are 

tenporary  intermediate  densities,  and  n,  8,  x,  X,  and  u  are  velocity- 
dependent  coefficients.  Here  x  and  X  are  diffusion  coefficients,  and  u  is 
the  antidiffusion  coefficient.  In  the  actual  algorithm,  <^+1  /0  is  corrected 

(hence  the  name  FCT)  to  a  value  $j+i/g  chosen  so  no  extrema  in  Pj  can  be 
enhanced  or  nev  ones  introduced  in  pj.  Previous  FCT  algorithms  had  0  *  C; 

the  widely  used  ETBFCT  and  related  algorithms  (Boris,  1976)  have  in  addition 
x  *  0.  If  ve  define  p,  to  be  sinusoidal  with  wave  number  k  on  a  mesh  with 
uniform  spacing  6x,  soJthat  p^  *  exp  (ij 6)  where  6  *  k6x,  then  the  nev 

density  array  satisfies 

Pj/p®  =  A  *  l-2i(n+8)sin8  +2(  x+X)  (cos  8-3  ) 

-  2u(cos  S-l)  [l-2insinS+2x(cr  -Z )  ] . 

From  A  we  can  determine  the  amplification  a  =  A  and  rel?  _.-e  phase  error  R  * 
(1/ eg)tan~1( -ImA/ReA)-l,  where  e  »  v6t/6x  is  the  Couran  umber. 

Expanding  in  powers  of  8  we  find 

a  «  1  +  a282  +  +  agS6  +  .  .  . 

R  ■  R0  +  R2s2  *  +  RgB6  +  .  .  . 

First-order  accuracy  entails  making  Rq  vanish,  which  requires  that  n  +  9  = 
e/2.  Second-order  accuracy  (og  *  0)  implies  that  y  =  k  +  X  -  e2/2. 
Analogously,  the  "reduced-phase-error"  property  Rg  ■  0  (Boris  and  Book, 

1976)  determines  u  *  (l-e2)/6,  thus  leaving  two  free  parameters.  One  of 
these  can  be  used  to  make  R^  vanish  also.  The  resulting  phase  error  R(8) 
is  small  not  only  as  M  0,  but  also  for  larger  values  of  8,  corresponding  to 
the  short  wavelengths  responsible  for  terraces  (Fig.  2) .  The  remaining 
parameter  n  can  be  chosen  to  relax  the  Courant  number  restriction  needed  to 
ensure  positivity  from  e  <  1/2  to  e  <  1.  When  coded,  these  changes 
necessitate  a  snail  increase  in  the  operation  count  of  ETBFCT  along  with  a 
small  increase. in  overhead  to  precalculate  the  two  new  arrays  of  velocity- 
dependent  transport  coefficients.  On  advection  tests,  the  new  algorithm 
completely  eliminated  terraces  (Fig.  3).  When  applied  to  the  coupled  systems 
of  gas  dynamic  equations,  it  produced  profiles  which  closely  approximate  the 
Riemann  solution  of  the  exploding  diaphragm  problem  (Fig.  U). 

The  second  approach  uses  a  rarefaction  flux  limiter  (RFL)  to  eliminate 
numerical  ripples  in  strong  rarefaction  waves.  This  approach  is  physically 
motivated.  Raw  anti-diffusive  fluxes  $j+]/2  are  limited  so  that 
the  slope  of  local  flow  field  profiles  decays  with  time  in  a  rarefaction 
wave.  In  effect,  additional  diffusion  is  left  in  the  field  to  maintain 
monotonicity  of  local  slopes.  For  mult i -material  calculations  a  "contact 
surface  sensor"  is  needed  to  detect  physical  discontinuities  and  shut  off  the 
RFL  locally. 

In  addition  we  found  that  some  care  was  required  when  applying 
generalized  continuity  equation  solvers  to  a  system  of  equations.  Truncation 
errors  of  the  various  equations  can  interact,  causing  undershoots  or 
overshoots  in  nonconvective  quantities  such  as  pressure.  We  found  that  it 
was  necessary  to  monotonize  derived  quantities  (pressure,  velocity)  before 
using  them  in  minimal-diffusion  transport  algorithms. 
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Hie  above  methodolog'  has  been  applied  be  a  series  of  best  problems 
inibiabed  by  a  spherical  high-explosive  (KE)  debonabior.  in  air.  Ar.  ideal 
Chapman -Jougueb  debonabion  was  used  bo  specify  bhe  inibial  condibions; 
afberbuming  was  neglecbed.  In  bhe  absence  of  reflecbing  surfaces, 
spherical  symmetry  is  maintained  and  bhe  calculabior.  remains  one- 
dimensional.  A  nonuniform  radial  grid  was  used  wibh  extremely  fine  zoning 
near  bhe  shock  fronb.  The  grid  was  moved  so  bhab  bhe  shock  remained 
approximabely  fixed  wibh  respeeb  to  the  mesh.  The  original  version  of  the 
FCT  algorithm  gave  rise  to  pronounced  terraces  in  bhe  rarefaction  region. 

This  would  have  rendered  any  two-dimensional  calculation  involving  shock 
diffraction  or  nonideal  effects  dubious.  The  techniques  described  here 
inproved  the  blast  wave  results  considerably.  The  decrease  in  phase  error 
reduced  terracing  dramatically. 

Next,  a  two-dimensional  (2D)  numerical  calculation  was  performed  bo 
simulate  one  of  Carpenter’s  (1974)  height -of -burst  experiments  which  used 
spherical  8-lb.  charges  of  PBX  9^04  at  51.6  cm.  The  previous  fine-zoned  ID 
calculation  was  used  to  initialize  the  problem.  It  was  mapped  onto  the  2D 
grid  Just  prior  to  the  onset  of  reflection.  The  solution  was  then  advanced 
in  time,  with  pressure  being  calculated  from  a  real-air  equation  of  state  and 
a  JVL  equation  of  state  for  the  combustion  products.  The  front  of  the  blasb 
wave  was  captured  in  a  finely  gridded  region  which  moved  outward  horizon¬ 
tally.  Special  care  was  taken  to  ensure  that  the  grid  moved  smoothly.  The 
resulting  solution,  particularly  the  curve  of  peak  overpressure  vs.  range, 
was  consistent  with  Carpenter's  experimental  data  (Fig.  8).  Although  this 
calculation  represents  a  reasonable  accurate  simulation  of  the  double-Mach- 
stem  region,  no  doubt  improvements  can  and  will  be  made  to  numerically  model 
such  phenomena. 
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Fig.  1  —  Rounded  half  circle  used  in  passive  scalar  advection  tests  (a)  initially, 
and  (b)  after  propagation  for  14  cycles  using  JPBFCT.  Note  that  terraces  form 
even,  as  here,  in  the  absence  of  comers  in  the  profile.  Tick  marks  indicate 
computational  zones  (N  ■  100). 


Fig.  2  —  Contour  plot  of  R(0,e)  for  new  mult  coefficient  FCT  algorithm 
Note  R  »  0  except  for  D  >  3  it  12.  The  relative  phase  error  vanishes  ex 
actiy  for  e  -  1/2  and  e  *  1 . 
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Fig.  3  —  (a)  Blowup  of  Fig.  1(a)  (dashed  line)  compared  with  (b)  same  profile  as 
computed  using  new  sixth  order-phase-accurate  FCT  algorithm.  Solid  traces  are 
exact  solutions. 


EXPLODING  DIAPHRAGM 
CYCLE  *120  TIME  *  0  4S2 


oo  : 

VELOCITY 

-3.85 

2500.0 

ENG  TOT 

25.0 

10 

DENSITY 


VELOCITY 

-3.97 

2500.0 

ENG  TOT 


Fig.  4 


25.0 


0  0  RADIUS  -  CM  -  100.0 


(b) 

(a)  Exact  and  (b)  computed  solution  of  exploding  diaphragm  problem 
(10-to-l  initial  density  jump,  100-to-l  initial  pressure  jump) 
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The  NRL  Flux-Corrected  Transport  code  FAST2D  is  used  to  mode)  the  conversion 
of  a  cylindrically  confined  surface  explosion  into  a  one-dimensional  blast  wave.  The  ideal 
gasdynamic  equations  are  used,  together  with  a  real-air  equation  of  state,  to  follow  the 
development  of  an  explosion  initialized  with  the  1-kton  standard  as  it  reflects  from  the 
cylindrical  boundary  and  propagates  upward.  Multiple  shocks  develop  and  the  solution 
quickly  (in  1-2  reflection  times)  approaches  the  one-dimensional  asymptotic  state 
described  by  the  Taylor-Sedov  similarity  solution. 
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CONVERSION  OF  A  CYLINDRICALLY  CONFINED  SURFACE  EXPLOSION 
INTO  A  ONE-DIMENSIONAL  BLAST  WAVE 


I.  INTRODUCTION 

In  connection  vith  the  Dense  Rack  (Close-Spaced  Rasing)  scheme,  it  has 
been  pointed  out  by  Latter  that  surface  environments  following  a  multiburst 
attack  could  be  considerably  more  difficult  to  survive  than  those  resulting 
from  single  bursts.  Detonation  simultaneously  over  all  the  shelters  of  a 
closely  spaced  array  surrounds  the  vicinity  of  an  individual  shelter  vith  high 
pressures,  allowing  the  explosion  to  vent  only  in  the  upward  direction. 

Hence,  at  least  in  the  interior  of  the  array,  extremely  high  pressures  are 
maintained  for  longer  times,  vith  the  result  that  total  impulses  applied  to 
ground  structures  are  nuch  greater  than  for  a  solitary  explosion. 

Latter  proposed  considering  the  following  idealized  situation:  identical 
explosions  are  initiated  simultaneously  in  an  infinite  regular  hexagonal  array 
at  (or  near)  a  uniform  level  plane  surface  (Fig.  1).  Construct  a  vertical 
plane  transverse  to  the  line  connecting  every  pair  of  neighboring  centers,  at 
the  midpoint  of  the  line.  By  symmetry,  the  six  planes  surrounding  an 
explosion  represent  perfectly  reflecting  surfaces  and  together  form  a 
hexagonal  parallelepiped  (a  cylinder  having  hexagonal  cross  section)  with  a 
vertical  axis.  Ibis  in  turn  can  be  approximated  as  a  right  circular  cylinder 
vith  the  same  axis  and  cross  sectional  area.  A  two-dimensional  r-z  code  can 
be  used  to  solve  this  idealized  problem  for  surface  and  air  bursts.  It  is 
clear,  however,  that  asymptotically  the  flow  becomes  one -dimensional,  evolving 
into  a  shock  tube  solution.  In  the  limit  where  this  asymptotic  state  has  been 
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reached  and  where  many  weapon  masses  of  air  have  been  swept  over,  the  one¬ 
dimensional  form  of  the  Sedov  point  blast  solution  should  apply.  The 
pressure  on  the  ground  as  a  function  of  time  t  is  given  by 

.  „  V  2/2  (1) 

t  ■ 
o 

where  p  is  the  air  density,  W  is  r,he  weapon  yield,  k  is  the  cylinder  area, 
o 

and  k  is  the  function  of  the  adiabatic  index  y  which  is  <  1.  This  is  tc  be 
compared  with  the  corresponding  result  for  the  case  of  a  spherical  free-field 
expansion,  which  is  of  the  form 


,  >  f  W 

p  *  k  Po  . 


As  may  be  seen,  the  pressure  given  by  Eq.  (l)  decreases  in  time  more 
slowly  than  that  of  Eq.  (2).  Although  the  idea  behind  this  so-called  "bomb- 
in-a-can"  model  is  simple  and  transparent,  it  leaves  several  important 
questions  unanswered.  First,  how  long  (how  many  shock  reflections)  does  it 
take  to  reach  the  asymptotic  state,  i.e.  ,  when  does  Eq.  (l)  become  a  good 
approximation?  Second,  to  what  extent  is  the  prediction  (l)  dependent  on 
idealizations  of  the  model,  e.g.  ,  perfect  symmetry,  perfect  simultaneity, 
neglect  of  scouring  and  entraining  of  dust,  etc?  finally,  what  are  the  long¬ 
term  tactical  consequences  in  this  model  of  a  multiburst  scenario— 
specifically,  what  can  be  said  about  the  shape  of  the  plume  produced,  how  nuch 
dust  is  raised  up  and  where  does  it  go,  what  sort  of  turbulence  is 


established  and  what  are  the  consequences  for  radar  transmission,  etc' 


In  this  report  we  describe  a  calculation  which  was  carried  out  to 
investigate  the  first  question,  regarding  the  approach  to  the  one -dimensional 
solution.  We  have  used  a  version  of  the  DNA  one  kiloton  (l-kton)  standard  to 


initialize  a  surface  burst  equivalent  to  18  scaled  megatons  (Mtor. 'i  ir.  a 
cylinder  of  radius  900  feet.  We  find  that  the  flov  is  close  to  or.e- 
dimensional  after  '*20  ms  (the  time  required  for  the  blast  to  reflect  from  t 
cylindrical  wall  and  return  to  the  origin,  considerably  less  than  might  hav 
been  expected. 

In  the  following  section  we  describe  in  detail  the  calculations  and  th 
results  obtained,  and  conclude  with  a  brief  discussion  of  what  we  intend  tc 


in  follow-on  work. 


2.  TWO  DIMENSIONAL  CALCULATION 
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The  development  of  the  explosion  can  be  viewed  in  terms  of  four  stages. 
These  are  (l)  the  spherical  free-field  expansion  prior  to  reflection  at  the 
outer  boundary;  (2)  multiple  distinct  reflections  off  both  boundaries,  during 
which  the  transition  to  one-dimensional  flow  takes  place;  (3)  approximately 
one  dimensional  expansion  after  reflecting  gas-dynamic  discontinuities  are  no 
longer  discernible;  (4)  a  gradual  transition  to  three-dimensional  flow  as  the 
finite  extent  of  the  array  comnunicates  itself  to  the  interior.  Stage  (l)  is 
not  contained  in  our  model;  since  the  weapon  mass  greatly  exceeds  the  mass  of 
air  contained  within  the  initial  explosion  radius,  the  nature  of  the  weapon 
and  its  constituents  ought  to  be  important  in  this  stage,  but  no  attempt  has 
been  made  to  improve  on  the  realism  of  the  1-kton  standard.  Likewise,  stage 
(4)  is  absent  from  the  model  (but  see  remarks  in  our  concluding  section 
below).  Instead,  we  attempted  to  model  stages  (2)  and  (3). 

The  calculation  was  performed  with  1Q0  zones,  each  10.2  cm  in  width  in 
the  radial  direction,  and  100  zones  in  the  axial  direction,  90  of  which  were 
20.4  cm  in  length  and  the  final  10  of  which  geometrically  increased  by 
increments  of  11 «  (Fig.  2a).  The  final  grid  size  was  10.2  m  by  22.4  m.  An 
equivalent  18  Mton  surface  burst  was  achieved  by  inserting  values  of  density, 
energy,  and  momentum  scaled  from  the  1  kton  standard.  Since  values  of  the  1 
kton  standard  were  only  available  from  10  meters  (an  18  Mton  surface  burst 
scales  to  about  eight  meters  after  we  double  the  yield  to  account  for  surface 
reflection),  we  created  initial  profiles  by  increasing  the  density,  internal 
energy,  and  momentum  by  a  factor  of  1.4.  A.s  a  result  the  initial  energy  was 
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3.1  x  1019  ergs  and  the  mass  vas  8. 2  x  106  g,  deposited  .into  a  hemisphere  of 
10  m  radius.  The  effect  of  this  scaling  vas  to  simulate  the  initiation  of  ar. 
event  with  the  appropriate  yield  and  weapon. mass,  but  with  profiles  (as  a 
function  of  distance  from  the  burst  site)  which  were  not  correct  in  detail. 
During  the  early  phases  of  the  problem  they  must  give  rise  to  discrepant 
behavior,  but  at  later  times  this  presumably  becomes  insignificant. 

The  numerical  simulation  vas  performed  using  the  JTRL  FAST2R  code  which 
utilizes  the  latest  version  of  Flux  Corrected  Transport  (FCT).  The  code  vas 
essentially  the  same  configuration  used  previously  for  lOU-ft  and  50-ft  1-kton 
HOB  calculations.  A  reflection  condition  was  inposed  on  the  left  and  bottom 
boundaries,  as  before,  and  on  the  right-hand  boundary  as  well.  In  addition  to 
the  diagnostics  employed  previously,  we  introduced  several  devised 
specifically  for  this  problem. 

Figure  2b  shows  plots  of  the  initial  profiles  of  the  mass  density, 
velocity,  pressure  and  energy  density.  Note  that  all  of  these  variables 
achieve  their  maxima  at  the  leading  shock.  Tb  convert  these  profiles  into 
their  counterparts  for  the  18  Mton  case,  it  is  necessary  to  dilate  by  a  factor 
of  (900  ft/10  ra)1/3  *  27. h.  Likewise,  the  time  scale  has  to  be  expanded  by 
the  same  factor. 

The  reflection  of  this  initial  blast  wave  from  the  outer  radius  of  the 
can  and  the  subsequent  evolution  of  the  system  are  depicted  in  Fig.  3.  This 
shows  pressure  contours  and  velocity  vectors  plotted  at  intervals  of  100 
cycles  (approximately  30-50  ys).  Although  the  time  intervals  vary  because  of 
the  changing  time  step,  these  plots  are  in  a  sense  equally  separated  in  terms 
of  the  "total  change"  in  the  system,  because  when  flow  and  sound  speeds  are 


large  the  timestep  decreases,  and  vice  versa.  Note  how  rapidly  the  reflected 

shock  propagates  to  the  left;  almost  invisible  at  cycle  201,  it  is  nearly 

halfway  across  the  mesh,  at  cycle  301.  It  travels  upward  almost  as  rapidly, 

reaching  an  altitude  of  10  m  by  cycle  601.  By  cycle  801  it  appears  to  have 

overtaken  the  primary  shock  wave,  which  is  considerably  distorted  after  cycle 

1501.  On  the  right  hand  boundary  the  reflection,  which  is  initially  regular, 

gradually  converts  to  Mach  reflection.  This  development  is  of  course 

analogous  to  the  transition  to  Mach  reflection  observed  in  the  BOB  case3.  It 

is  hard  to  see  exactly  when  transition  occurs,  but  a  Mach  stem  is  clearly 

visible  in  the  pressure  plots  at  cycle  (1*0)  and  thereafter.  If  we  examine 

Carpenter's  results  regarding  the  boundary  between  regular  and  Mach  reflect io 
in  real  air  at  Mach  numbers  M  >  10,  we  see  that  the  angle  at  which  transition 

occurs  depends  on  M  only  weakly  and  is  approximately  equal  to  15°.  The  Mach 
stem  ought  thus  to  have  first  formed  around  cycle  701  to  801,  but  the  region 
of  reflection  is  net  sufficiently  well  resolved  to  show  it. 

In  Fig.  1*  another  view  of  the  same  time  sequence  is  shown.  These  surfac 
pressure  plots  are  particularly  well  suited  for  showing  shock  and  rarefaction 
waves  propagating  through  the  sytem.  Since  the  pressures  are  interpolated 
onto  a  regular  mesh  before  plotting,  peak  values  are  lowered  by  as  much  as  a 
factor  of  two.  This  reduces  contrast  somewhat  and  makes  some  features  appear 
to  move  up  and  down  erratically  in  time,  but  the  overall  morphology  is  very 
clear.  We  see,  for  example,  how  the  shock  reflected  radially  from  the  origin 
appears  to  run  out  of  gas  at  about  cycle  901,  while  a  train  of  waves 
reverberates  down  from  the  top.  One  of  these  finally  reaches  the  radially- 
propagated  ground  shock  at  about  cycle  1801  and  begins  to  push  it  along. 
Throughout  this  time  the  pressures  in  the  lower  right  corner  of  the  system 
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remain  lov.  Meanwhile,  the  shock  front  rising  toward  the  top  cf  the  system, 
which  is  madly  churning  around  at  cycle  1201,  becomes  smoother  as  shocks 
propagate  back  down  from  it. 

Figure  5  shows  the  pressure  recorded  at  a  series  of  sensors  located 
across  the  bottom  of  the  can,  ranging  from  the  burst  center  (a'  to  the  outer 
edge  (1).  Note  that,  although  the  initial  profiles  have  a  sharp  (one-cone) 
discontinuity  at  the  leading  edge,  the  pressure  peak  at  sensor  £  takes  about 
TO  us  to  build  up  to  the  maximum.  This  is  about  a  hundred  times  as  long  as  it 
takes  the  shock  to  propagate  across  one  zone.  It  is,  however,  net  out  of 
line  with  the  time  required  for  the  air  behind  the  shock  (velocity  1.5  x  1C6 
cm/s)  to  cross  a  zone,  about  7  us.  The  density  (and  other  fluid  variables) 
can  build  up  to  their  post-reflection  values  when  the  air  from  ( y+1) / ( y-1) 

*  10  cells  has  been  compressed  into  the  cell  closest  to  the  wall. 

Using  Carpenter’s  **  reflection  factors  for  normally  incident  shocks,  we 
calculate  that  the  free-field  pressure,  initially  3.5  kbar,  should  reflect  up 
to  about  1*0  kbar,  which  is  to  be  cotrpared  with  the  value  of  25  kbar  we 
actually  obtained.  If  we  look  at  sensors  closer  to  the  burst,  site,  we  see 
that  the  pressure  peaks  get  continuously  sharper,  finally  turning  into  an 
almost  discontinuous  spike  at  the  origin,  as  the  imploding  shock  attempts  to 
develop  into  a  singularity. 

The  information  in  these  pressure  histories  is  capsulized  in  Fig.  6, 
which  shows  the  maximum  pressure  recorded  during  the  course  of  the  calculation 
as  a  function  of  sensor  location.  Note  that  the  peak  at  r  *  0  is  nearly  as 
high  as  that  at  the  periphery,  in  contrast  with  results  reported  by  TVatt5. 
This  may  be  attributable  to  the  lower  resolution  he  used,  or  to  the  greater 
numerical  dissipation  in  the  HULL  code.  Only  one  peak  was  detected  at  each 
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sensor  location.  Our  calculation  was  run  for  less  than  1  ms,  or  (scaled  to  IP 
Mton)  rather  less  than  the  30  ms  Fyatt  required  to  see  a  second  peak  at  the 
periphery  when  the  re-reflected  shock  impinges  there ;  the  two  calculations 
are  thus  not  inconsistent  in  this  regard. 

The  approach  to  one-dimensional  motion  may  be  inferred  from  the  contour 
plots  of  Fig.  U,  where  it  is  seen  that  the  leading  edge  of  the  upward- 
propagating  blast  wave  becomes  more  and  more  horizontal,  and  from  the  velocity 
vector  plots  of  Fig.  3,  which  show  that  the  flow  becomes  increasingly  vertical 
and  that  velocities  well  behind  the  upper  front  tend  to  die  away.  A  more 
quantitative  comparison  is  possible  if  we  compute  the  average  pressure  on  the 
bottom,  R 


2ir  /  p(r,0)rdr 


n 


as  a  function  of  time.  The  result  is  shown  as  the  solid  line  in  Fig.  7.  For 
comparison,  formula  (l)  has  been  added  (broken  line),  using  yield  W  *  1.1  kton 
and  area  A  *  lOOir  m2,  with  y  =  1.2.  The  agreement  between  the  two  is 
striking,  particularly  since  we  expect  to  observe  it  only  asymptotically .  It 
is  clear  that  the  one-dimensional  approximation  becomes  valid  very  early, 
probably  because  the  reverberating  shocks  behind  the  leading  front  of  the 
blast  propagate  so  rapidly  in  the  hot  medium. 


o 
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3.  CONCLUSIONS 

We  have  successfully  modelled  the  "bomb-in-a  car.”  using  FAST2L,  a  general- 
purpose  code  which  required  no  special  modification  for  this  problem. 

Although  we  only  ran  the  calculation  out  to  an  altitude  of  about  20  r.  (twice 
the  system  radius),  the  sliding  grid  (continuous  adaptive  recone)  capability 
of  FAST2D  allows  us  to  continue  the  run  as  long  as  desired  by  stretching  the 
mesh  in  the  vertical  direction. 

Our  results  show  that  the  asymptotic  one-dimensional  state  is  approached 
very  rapidly,  apparently  because  the  multiply  reflected  shocks  propagate  much 
faster  than  the  original  blast  wave.  Two  or  three  reflections  across  the 
radius  of  the  system  effectively  equilibrate  the  flow  and  relax  it  to  a  state 
of  vertical  expansion.  The  simulation  reveals  a  complex  pattern  of 
reverberations,  with  reflection  occurring  off  all  the  boundaries  and  regions 
where  conditions  are  highly  rronuniform. 

A  number  of  modifications  are  required  to  make  these  results  applicable  to 
an  actual  tactical  situation.  Perhaps  the  most  important  is  the  inclusion  of 
dust  scoured  up  from  the  ground  and  entrained  in  the  wind  fields  following  the 
blast*  A  realistic  description  would  require  not  only  that  the  dust  mass  load 
the  air,  but  that  air  and  dust  -be  permitted  to  exchange  momentum  and  energy. 
The  bottom  boundary  conditions  should  also  be  changed  to  model  the  energy  that 
goes  into  cratering  and  scouring. 

If  the  calculation  is  to  be  continued  to  late  times,  atmospheric 
stratification  should  be  included.  The  effects  of  water  vapor  and  turbulence 
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should  also  he  modelled.  Also,  venting  of  the  explosion  in  the  horizontal 
direction  (in  the  case  of  an  RV  exploding  near  the  edge  of  a  finite-sized 
array  of  shelters,  or  for  an  explosion  in  the  interior  of  the  array  which  has 
begun  to  feel  the  pressure  in  the  outer  explosions  drop  as  they  vent  outward) 
can  be  modelled  by  having  the  radius  R  increase  as  a  function  of  height.  'Of 
course,  airblast  situations  in  which  the  explosion  originates  above  the  ground 
can  also  be  investigated  (in  which  case  Mach  stems  could  appear  on  the  ground 
as  well  as  the  sides  of  the  system). 

Although  some  of  these  extensions  require  nontrivial  code  modifications, 
there  is  no  reason  in  principal  why  the  present  results  could  not  be  augmented 
and  refined  greatly  by  additional  calculations. 


Fig.  1  —  Schematic  of  unbounded  hexagonal  array  of  shelters  with  1800  ft 
separations  on  a  perfectly  reflecting  plane.  Shown  is  the  cross  section  of  the 
hexagonal  paxallelpiped  formed  by  the  planes  bisecting  the  lines  connecting 
a  given  shelter  with  its  six  nearest  neighbors  (reflection  planes),  along  with 
the  coaxial  cylinder  having  the  same  cross  sectional  area. 


PRESSURE  CONTOURS  VELDC I TY  VECTORS 


Pressure  contours  and  velocity  vectors  calculated  using  FAST2I),  shown  at  intervals  of  100  cycles 
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Fig.  3b  —  Pressure  contours  and  velocity  vectors  calculated  using  FAST2D,  shown  at  intervals  of  100  cycles 


* 

o 

« 

a.  •  o 

Ul  UIO 
UJ  J  gj 
z  u  w 

UJ  >- 

o  u 


£  L 

OL  •  o 
uj  u  a 
UJ 

z  u 

UJ  >- 

a  u 


wv 


MAX  VELDC I TY  '5.21*10 

Fig.  3d  —  Pressure  contours  and  velocity  vectors  calculated  using  FAST2D,  shown  at  intervals  of  1 00  cycles 


MAX  VELOCITY  •  i. 56*  10 

Pressure  contours  and  velocity  vectors  calculated  using  FAST2D,  shown  at  intervals  of  100  cycles 
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Fig.  3h  —  Pressure  contours  and  velocity  vectors  calculated  using  FAST2D,  shown  at  intervals  of  100  cycles 
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Fig.  3j  —  Pressure  contours  and  velocity  vectors  calculated  using  FAST2D,  shown  at  intervals  of  100  cycles 
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Fig.  3k  -  Pressure  ran  tours  and  velocity  vectors  calculated  using  FAST2D,  shown  at  intervals  of  100  cycles 
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Fig.  3p  -  Pressure  contours  and  velocity  vectors  calculated  using  FAST2D,  shown  at  intervals  of  100  cycles 
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Fig.  4b  —  Pressure  contours  of  Fig.  3,  shown  in  orthographic  projection. 
The  pressure  scale  here  is  22  kbar. 


Fig.  4c  —  Pressure  contours  of  Fig.  3,  shown  in  orthographic  projection 
The  pressure  scale  here  is  22  kbar. 


Fig.  4d  -  Pressure  contours  of  Fig.  3,  shown  in  orthographic  projection. 
The  pressure  scale  here  is  22  kbar. 


Fig.  4h  —  Pressure  contours  of  Fig.  3,  shown  in  orthographic  projection. 
The  pressure  scale  here  is  22  kbar. 
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Pig.  4i  —  Pressure  contours  of  Fig.  3,  shown  in  orthographic  projection 
The  pressure  scale  here  is  22  kbar. 
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Pig.  4j  —  Pressure  contours  of  Fig.  3,  shown  in  orthographic  projection 
The  pressure  scale  here  is  22  kbar. 
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Fig.  4q  —  Pressure  contours  of  Fig.  3,  shown  in  orthographic  projection 
The  pressure  scale  here  is  22  kbar. 


Fig.  5a  —  Pressure  histories  as  measured  by  sensors  located  at  radii  of  0.0m 


Fig.  5b  Pressure  histories  as  measured  by  sensors  located  at  radii  of  0.035m 
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Fig.  5c  —  Pressure  histories  as  measured  by  sensors  located  at  radii  of  0.075m 
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Fig.  5d  —  Pressure  histories  as  measured  by  sensors  located  at  radii  of  1 ,25m 
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Fig.  5e  —  Pressure  histories  as  measured  by  sensors  located  at  radii  of  2.0m 
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Fig.  5f  —  Pressure  histories  as  measured  by  sensors  located  at  radii  of  3.0m 
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Fig.  5g  —  Pressure  histories  as  measured  by  sensors  located  at  radii  of  4 ,5m 
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Fig.  5h  —  Pressure  histories  as  measured  by  sensors  located  at  radii  of  6.05m 
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Fig.  6  —  Maximum  recorded  station  pressure  in  kbar  as  function 
of  station  radial  location 
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The  flow  fields  that  result  from  nuclear  and  high  explosive  (HE)  detonations  are  qualitatively  alike 
but  quantitatively  different.  Consequently,  care  must  be  exercised  in  carrying  over  conclusions  drawn 
from  measurements  of  HE  tests  to  nuclear  explosions.  The  usefulness  of  HE  explosions  for  simulating 
nuclear  airblast  is  predicated  on  the  fact  that  after  reaching  5-6  times  the  initial  radius  the  flow  field  looks 
like  that  produced  by  a  point  source  and  produces  shock  overpressures  simular  to  those  in  the  nuclear 
case.  Numerical  simulations  of  airblast  phenomena  have  been  carried  out  using  one-  and  two-fluid 
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Flux-Corrected  Transport  hydrody codes  in  one  and  two  dimensions.  The  principal  difference  in  the 
free-field  solutions  are  the  presence  in  the  HE  case  of  a  contact  discontinuity  between  air  and  HE 
products  and  of  a  backward-facing  shock  behind  it.  Temperatures  in  the  nuclear  fireball  are  initially 
three  orders  of  magnitude  higher;  correspondingly ,  the  density  minimum  at  the  center  of  the  fireball  is 
much  broader  and  deeper.  When  the  blast  wave  in  a  nuciear  heigh t-of-burst  (HOB)  situation  undergoes 
regular  reflections  from  the  ground  only  one  peak  develops  in  the  overpressure,  and  the  reflected  wave 
propagates  upward  rapidly  through  the  hot  underdense  fireball.  In  the  HE  case  part  of  the  upward-moving 
reflected  wave  is  reflected  downward  at  the  contact  surface,  producing  a  second  pressure  peak  on  the 
ground,  while  the  shock  transmitted  through  the  contact  surface  propagates  slowly  upward.  After  tran¬ 
sition  to  Mach  reflection  other  differences  appear.  At  late  times  following  shock  breakaway  the  nuclear 
fireball,  unlike  the  HE  fireball,  appears  to  develop  a  Rayleigh  -Taylor  instability  along  its  lower  edge  below 
the  HOB.  The  vortices  (both  forward  and  reverse)  are  stronger  and  form  earlier.  This  has  important 
consequences  for  fireball  rise  and  for  dust  entrainment  and  transport  to  high  altitudes. 
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The  flow  fields  that  result  from  nuclear  and  high 
explosive  (HE)  detonations  are  qualitatively  alike  but 
quantitively  different.  Consequently,  care  must  be  exercised 
in  carrying  over  conclusions  drawn  from  measurements  of  HE 
tests  to  nuclear  explosions.  The  usefulness  of  HE  explosions 
for  simulating  nuclear  airblast  is  predicated  on  the  fact 
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field  looks  like  that  produced  by  a  point  source  and  produces 
shock  overpressures  similar  to  those  in  the  nuclear  case. 
Numerical  simulations  of  airblast  phenomena  have  been  carried 
out  using  one-  and  two-fluid  Flux-Corrected  Transport  hydro¬ 
codes  in  one  and  two  dimensions.  The  principal  differences 
in  the  free-field  solutions  are  the  presence  in  the  HE  case 
of  a  contact  discontinuity  between  air  and  HE  products  and  of 
a  backward-facing  shock  behind  it.  Temperatures  in  the 
nuclear  fireball  are  initially  three  orders  of  magnitude 
higher;  correspondingly,  the  density  minimum  at  the  center  of 
the  fireball  is  much  broader  and  deeper.  When  the  blast  wave 
in  a  nuclear  height-of-burst  (HOB)  situation  undergoes 
regular  reflection  from  the  ground  only  one  peak  develops  in 
the  overpressure,  and  the  reflected  wave  propagates  upward 
rapidly  through  the  hot  underdense  fireball.  In  the  HE  case 
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part  of  the  upward-moving  reflected  wave  is  reflected 
downward  at  the.  contact  surface,  producing  a  second  pressure 
peak  on  the  ground,  while  the  shock  transmitted  through  the 
contact  surface  propagates  slowly  upward.  After  transition 
to  Mach  reflection  other  differences  appear.  At  late  times 
following  shock  breakaway  the  nuclear  fireball,  unlike  the  HE 
fireball,  appears  to  develop  a  Ray leigh-Taylor  instability 
along  its  lower  edge  below  the  HOB.  The  vortices  (both 
forward  and  reverse)  are  stronger  and  form  earlier.  This  has 
important  consequences  for  fireball  rise  and  for  dust 
entrainment  and  transport  to  high  altitudes. 


Section  1 
INTRODUCTION 

In  this  paper  we  describe  a  series  of  calculations 
carried  out  as  part  of  the  ongoing  NRL  effort  aimed  at  study¬ 
ing  airblast  effects.  The  phenomena  of  chief  interest  to  us 
include  the  following:  peak  overpressures  and  pressure 

histories  on  the  ground  as  functions  of  yield,  range,  and 
height  of  burst  (HOB),  both  at  early  times  (prior  to  and 
during  transition  to  Mach  reflection)  and  at  late  times 
(after  shock  breakaway,  with  peak  pressure  in  the  range  of 
tens  of  psi );  velocity  fields,  particularly  those  associated 
with  the  toruses  (both  forward  and  reverse)  in  the  neighbor¬ 
hood  of  the  rising  fireball;  and  the  distribution  of  dust 
lifted  off  the  ground  by  the  winds  and  the  structure  of  the 
cloud  at  the  time  of  stabilization.  We  are  interested  in 
comparing  the  nuclear  and  HE  cases,  and  learning  how  much 
they  differ  from  one  another.  Our  motivation  is  to  determine 
the  extent  to  which  HE  tests  can  simulate  events  in  a  nuclear 
heigh t-of-burst  situation. 
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The  technique  we  have  employed  for  this  purpose  is 
numerical  modeling.  One-  and  two-fluid  hydrocodes  based  on 
the  Flux-Corrected  Transport  { FCT )  shock-capturing  techni¬ 
ques  A  have  been  used  to  simulate  airblast  phenomena  in  one 
and  two  dimensions.  FCT  refers  to  a  class  of  state-of-the- 
art  fluid  computational  algorithms  developed  at  NRL  in  the 
course  of  the  past  ten  years  with  supersonic  gas-dynamic 
applications  expressly  in  mind.  Simply  put,  our  procedure  is 
to  model  a  one-kton  nuclear  burst  and  its  600-ton  chemical 
equivalent,  both  at  a  HOB  of  50  m,  and  compare  the  results. 
In  order  to  validate,  initialize,  and  interpret  these  2D 
simulations,  a  number  of  ancillary  calculations  (mostly  ID) 
were  undertaken.  The  results  are  most  conveniently  exhibited 
in  terms  of  plots  of  peak  overpressure  vs  range  and  time, 
station  histories,  contour  plots  of  combustion  product  and 
total  density,  velocity  vector  plots,  and  tracer  particle 
trajectories.  Examples  of  these  are  presented  to  illustrate 
our  results  and  conclusions. 
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The  plan  of  the  paper  is  as  follows:  In  the  next 
section  we  discuss  our  numerical  techniques  and  validation 
procedures.  In  Section  3  we  discuss  the  free-field  (ID) 
solution  and  indicate  the  salient  differences  between  nuclear 
and  HE  cases.  Section  4  describes  the  2D  HOB  calculations 
done  for  the  HE  and  nuclear  cases.  In  Section  5  we  summarize 
our  conclusions  and  discuss  their  domain  of  validity.  We 
find  that  simulation  of  nuclear  explosions  by  HE  has  distinct 
limitations,  particularly  at  early  times,  in  the  fireball  and 
transition  regions,  and  in  the  details  of  the  dust  scouring 
process. 
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Section  2 

NUMERICAL  TREATMENT 


FCT  is  a  finite-difference  technique  for  solving 
the  fluid  equations  in  problems  where  sharp  discontinuities 
arise  (e.g.,  shocks,  slip  surfaces  and  contact  surfaces).1 
It  modifies  the  linear  properties  of  a  second-  (or  higher-) 
order  algorithm  by  adding  a  diffusion  term  during  convective 
transport,  and  then  subtracting  it  out  "almost  everywhere"  in 
the  antidiffusion  phase  of  each  time  step.  The  residual 
diffusion  is  just  large  enough  to  prevent  dispersive  ripples 
from  arising  at  the  discontinuity,  thus  ensuring  that  all 
conserved  quantities  remain  positive.  FCT  captures  shocks 
accurately  over  a  wide  range  of  parameters.  No  information 
about  the  number  or  nature  of  the  surfaces  of  discontinuity 
need  be  provided  prior  to  initiating  the  calculation. 

The  FCT  routine  used  in  the  present  calculations, 
called  JPSFCT  (an  advanced  version  of  ETBFCT)2,  consists  of  a 
flexible,  general  transport  module  which  solves  1-D  fluid 
equations  in  Cartesian,  cylindrical,  or  spherical  geometry. 
It  provides  a  finite-difference  approximation  to  conservation 
laws  in  the  general  form: 

tt  /  >dV  -  -  /_  o(u-u)-dA  +  f  xdA,  (1) 

5V(  t )  5A(  tT  ~  SA(  t ) 

where  s  represents  the  mass,  momentum,  energy  or  mass  species 
in  cell  5V(t),  u  and  Ug  represent  the  fluid  and  grid  velo¬ 
cities,  respectively,  and  t  represents  the  pressure/work 
terms.  This  formulation  allows  the  grid  to  slide  with 
respect  to  the  fluid  without  introducing  any  additional 
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numerical  diffusion.  Thus,  knowing  w.nere  the  features  of 
greatest  interest  are  located,  one  can  concentrate  fine  zones 
where  they  will  resolve  these  features  most  effectively  as 
the  system  evolves. 

The  same  transport  routine  is  employed  in  the  2D 
r-z  code  (called  FAST 2D )  via  coordinate  splitting.  A  Jones- 
Ailkins-L.ee  ( JV»  L )  equation  of  state  (EOS)  was  used  for  the 
detonation  products  and  a  real-air  EOS  was  used  outside  the 
HE-air  interface.3  The  routine  was  written  in  the  form  of  a 
table  lookup,  using  interpolation  with  logarithms  to  the  base 
16  computed  by  means  of  logical  shifts.1*  By  thus  taking 
account  of  the  architecture  of  the  machine  (in  these  calcula¬ 
tions,  a  32-bit-word  two-pipe  Texas  Instruments  ASC)  it  was 
possible  to  generate  very  efficient  vector  code,  decreasing 
the  time  required  for  EOS  calculations  to  a  small  fraction  of 
that  required  for  the  hydro.  The  EOS  specifies  pressure  as  a 
function  of  density  and  internal  energy.  In  mixed  cells  the 
combined  pressure  was  calculated  according  to  Dalton's  law. 


For  the  HE  calculations  the  initial  conditions  were 
taken  to  be  the  self-similar  flow  field  corresponding  to  a 
spherical  Chapman-Jouguet  detonation  at  the  time  the  detona¬ 
tion  wave  reaches  the  charge  radius  (Fig.  iH.  This  was 
propagated  with  the  ID  spherical  code  until  the  detonation 
front  attained  a  radius  just  smaller  than  the  HOB,  at  which 
time  the  solution  was  laid  down  on  the  2D  mesh  (Fig.  ,2).  .  The 
nuclear  calculation  was  initialized  with  the  1-kton  standard5 
with  the  same  initial  radius. 

The  boundary  conditions  were  chosen  to  enforce 
perfect  reflection  on  the  ground  and  on  the  axis  of  symmetry 
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[  ( a  $  / .  ).  =  0,  where  :  =  p,o,v  ,  ana  v.  =  01,  where  "t  ana 

an  oc  dc 


’n"  denote  tangential  and  normal  components,  respectively, 


and  outflow  on  the  outer  and  the  too  boundaries  ! (do/dn),  = 

t  n  i  -be 

C,  where  5  *  p  ,p ,  v,  v  j . 


For  the  2D  calculations  the  mesh  was  typically 
~100  x  100  .  Fixed  gridding  was  used  to  minimize  numerical 
errors.  These  zones  were  2.1  m  x  2.1  m.  For  the  late  time 
calculations,  a  fixed  mesh  with  100  zones  in  the  radial  and 
200  zones  in  the  vertical  direction  was  used,  with  all  cells 


of  dimension  4.2  m  x  4.2  m. 


Section  3 


FREE-FIELD  SOLUTION 


The  well-known  Sedov  similarity  solution6  for  a 
point  blast  consists  of  a  strong  shock  (post-shock  pressure 
much  larger  than  ambient  pressure)  followed  by  a  rarefaction 
wave  (Fig.  3).  The  density  distribution  is  extremely  concave 
and  approaches  zero  at  the  origin.  The  pressure  approaches  a 
constant  as  r  +  0,  so  that  the  temperature  diverges  strongly. 
The  profiles  of  the  1-kton  standard  solution  (Fig.  4)  are 
qualitatively  similar,  the  temperature  being  essentially  -  flat 
in  the  fireball  region. 


The  solution  used  to  initialize  the  HE  problem, 
however,  contains  a  number  of  features  which  are  absent  in 
the  other  solutions.  These  are  most  conspicuous  in  the 
density  profile  (Fig.  2),  which  exhibits  a  contact  discon¬ 
tinuity  between  the  HE  products  and  shocked  air,  a  secondary 
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shock  facing  inward  within  the  detonation  products,  and  a 
gentle  maximum  near  the  origin.  Because  of  this  last 
feature,  temperatures  are  about  three  orders  of  magnitude 
smaller  in  the  fireball  than  in  the  nuclear  case,  and  are 
nowhere  divergent. 


It  follows  that  the  speed  of  sound  in  the  nuclear 
fireball  is  much  greater  than  in  the  HE  fireball.  This  has 
two  immediate  consequences,  one  physical  and  one  numerical. 
The  first  is  that  shocks  propagating  through  the  nuclear 
fireball  travel  much  faster.  The  second  is  that  the  upper 
limit  on  the  computational  timestep,  set  by  the  Courant 


criterion 


max 


vj+c 


.  )  <  1 


which  usually  is  determined  by  conditions  in  the  fireball,  is 

much  smaller  relative  to  the.  shock  time  scale  t  ■ 

s 

HOB/v  ,  ,  in  the  nuclear  case  than  in  the  HE  case. 

snocK 

As  a  result,  even  though  the  leading  shock  is  well 
resolved  (over  -  2  zones)  in  the  free-field  solution,  the 
process  of  reflection  even  at  ground  zero  (where  the  shock  is 
incident  normally)  takes  hundreds  of  timesteps.  Coupled  with 
the  property  of  FCT  (known  as  "clipping")  which  makes  the 
points  of  all  sharply-peaked  profiles  tend  to  flatten  out 
until  they  are  >  3  zones  across,  this  makes  the  rise  time  of 

the  reflected  shock  quantities  ~  10  times  longer  than  that  of 

the  incident  shock.  This  can  be  seen  using  a  ID  spherical 
model  calculation  (Fig.  4),  as  well  as  in  oblique  reflection 
in  2D.  This  spreading  is  a  problem  only  while  the  shock  is 
in  the  immediate  vicinity  of  the  reflecting  boundary.  After 


the  reflected  shock  has  propagated  a  few  zones  back  into  the 
interior  of  the  mesh,  the  profiles  steepen  up  and  assume 
their  correct  forms  (this  has  been  shown  by  rerunning  the 

calculation  with  a  refined  mesh  and  comparing  the  peak  values 
after  reflection  with  those  predicted  theoretically). 

When  the  reflected  shock  begins  propagating  back  to 
the  origin  it  encounters  drastically  different  conditions  in 
the  HE  and  nuclear  cases.  In  the  former,  it  strikes  the 

contact  discontinuity  where  it  is  partly  transmitted  and 
partly  reflected.  The  reflected  shock  then  proceeds  outward 
until  it  reaches  the  ground  (the  end  of  the  grid  in  the  ID 
calculation),  producing  the  second  peak  in  the  station 
history  shown  in  Fig.  5.  In  contrast,  the  shock  wave 
reflected  from  the  ground  passes  unhindered  through  the 

fireball  at  high  speed  until  it  reaches  the  upper  boundary  of 
the  fireball  whereupon  it  reflects  back. 

As  we  shall  show,  it  is  primarily  through  these 
reverberating  shock  waves  and  the  wind  pattern  they  set  up 

that  the  HE  and  nuclear  HOB  airblasts  differ. 


Section  4 

2D  SIMULATION  OF  AIRBLAST 


The  yield  and  HOB  (  600  tons  and  166  ft,  respec¬ 
tively)  were  chosen  to  equal  the  values  used  in  the  Direct 
Course  experiment,  which  we  are  simulating.  The  Chapman- 
Jouguet  parameters  used  to  initialize  the  spherical  free- 
field  calculation  were  taken  to  be  those  for  the  NH4N03-fuel 
oil  (ANFO)  mixture  used  as  the  explosive. 


Figures  5(a)-(c)  show  the  contours  of  HE  density 
and  internal  energy  per  unit  mass  and  the  velocity  arrow  plot 
at  t*Q,  just  before  the  reflection  at  ground  zero  occurs. 
Figures  6 ( d } - ( f )  show  the  corresponding  plots  54  ms  later, 
while  Figs.  6(g)— (i)  show  them  after  245  ms.  Note  the 
reflected  shock  proceeding  upward,  reflecting  again  off  the 
fireball,  and  propagating  back  in  a  downward  and  outward 
direction.  The  interaction  of  this  shock  with  the  radially 
inward  flow  near  -.the  ground  generates  the  reverse  vortex, 
which  is  clearly  seen  in  Fig.  6(i).  Note  also  the  positive 
vortex  forming  near  the  top  of  the  grid  in  the  same  plot. 
The  latter  results  when  the  upward-propagating  reflected 
shock  interacts  with  the  radially  outward  flow  near  the  top 
of  the  fireball;  it  is  not  produced  by  the  buoyant  rise  of 
the  fireball,  which  at  these  early  times  has  scarcely  begun. 

To  look  at  the  evolution  of  the  fireball  at  late 
times,  we  reinitialized  on  a  larger,  coarser  grid,  represent¬ 
ing  a  cylinder  400  m  in  radius  and  300  m  high.  The  first  300 
cycles  approximately  reproduce  the  early-time  results.  The 
spherical  shock  breaks  away  and  leaves  the  mesh.  The  flows 
remaining  on  the  grid  are  now  subsonic  everywhere.  Then  the 
fireball  begins  to  rise  and  the  subsequent  development  is  due 
to  the  combination  of  buoyant  rise  and  the  action  of  the 
vortices  set  up  by  the  early  shocks. 

Figures  7(a)-(b)  show  the  reaction  product  density 
and  velocities  at  0..93  s.  Note  the  "toe"  reaching  out  along 
the  ground  and  the  bulge  near  the  bottom  of  the  HE  product 
density  produced  by  the  constructive  interference  of  forward 
and  reverse  vortices.  These  features  are  accentuated  with 
the  passage  of  time;  in  Figs.  7(c) -(d)  (t  =  2.70  s),  they  are 
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Figure  6f 
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ever,  clearer.  The  ciouc  has  become  quire  elongated  verti¬ 
cally  and  shows  a  distinct  mushroom  shape.  Development  slows 
as  the  fireball  cools  and  velocities  diminish.  By  7.34  s 
(after  2600  timesteps )  the  cloud  is  almost  at  600  m.  Figures 
7(e)-(f)  show  its  form  at  this  time.  Note  that  the  maximum 
velocity  is  now  145  m/s. 


Figure  8a,  which  shows  the  trajectories  of 
passively  advectea  tracer  particles  over  the  time  interval 
1.8  sec  to  3.97  sec,  displays  the  vortices  very  clearly. 
Figure  8b  shows  the  particle  paths  for  the  time  interval  3^97 
sec  to  7.34  sec.  Notice  that  there  are  four  vortices  visible 
in  the  plot:  two  positive  and  two  reversed.  The  additional 
small  vortices  are  apparently  a  consequence  of  entrainment  by 
the  major  ones.  As  far  as  we  know,  their  existence  has  not 
been  noted  previously. 

When  we  repeat  the  calculation  with  nuclear  initial 
conditions,  several  differences  appear  at  a  very  early  stage. 
The  reflected  shock  propagates  upward  rapidly  through  the 
much  hotter  fireball  and  reverberates  more.  The  maximum  flow 
speeds  (as  opposed  to  sound  speeds)  are  smaller,  a  difference 
which  persists  to  late  times.  Although  the  shock  radius  as  a 
function  of  time  is  essentially  the  same,  the  rarefaction 
wave  moves  faster  as  the  deeper  density  well  gets  filled  in. 

Figures  9(a)-(c)  show  plots  analogous  to  those  of 
Figures  6(g)  —  (i) ;  by  this  time  it  is  clear  that  much  of  the 
early  difference  in  the  density  profiles  is  washed  out  as 
pressure  begins  to  relax  to  ambient.  (The  pressure  differs 
from  ambient  by  <  5%  everywhere,  so  that  pressure  contour 
plots  are  not  very  informative.)  Note,  however,  the  indenta¬ 
tions  that  appear  on  the  underside  of  the  internal  energy 
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contours  at  trie  base  of  tne  sterr,  [Fig.  9(b)j.  Tnese  are 
absent  in  the  corresponding  HE  plot,  Fig.  o(g).  We  conjec¬ 
ture  that  they  are  the  signature  of  a  fluid  instability, 
possible  Ravleigh-Tavlor .  The  idea  is  that  the  air  sucked  in 
by  the  (forward)  vortex  at  the  bottom  of  the  fireball  is  much 
denser  than  the  fireball  itself.  In  running  into  the  latter 
it  sets  up  the  classic  condition  for  Ravleigh-Taylor  insta¬ 
bility  (direction  of  effective  gravity  and  density  gradient 
are  opposed ) . 

It  is  clear  that  the  major  qualitative  differences 
between  the  HE  and  nuclear  cases  persist  longest  in  the 
velocity  plots.  This  is  not  surprising,  as  the  circulation 
patterns  represented  by  the  vortices  have  essentially 
infinite  lifetimes  in  the  absence  of  viscosity.  we  have  run 
both  -  nuclear  and  HE  cases  out  to  stabilization  (not  shown 
here)  and  have  shown  that  there  are  qualitative  differences 
in  velocity  plots  to  the  very  end.  At  all  times  t>0  the  peak 
flow  velocity  in  the  HE  case  exceeds  that  in  the  comparable 
nuclear,  result.  This  is  a  reflection  of  the  fact  that  the 

Chapman-Jouguet  solution  at  a  radius  of  10  m  has  a  pressure 
peak  of  52  kbar,  vs  3  kbar  for  the  1-kton  standard  at  the 
same  radius.  The  means  that  the  former  starts  out  with  much 

more  violent  motion,  i.e.,  fluid  velocities  an  order  of 

magnitude  larger.  In  point  of  fact,  the  HE  case  does  not 

closely  resemble  a  point  source.  At  initialization  the 
nuclear  profiles  have  ~  6%  of  the  yield  in  kinetic  energy. 
This  fraction  increases  to  a  maximum  of  ~  15%,  then 

decreases.  In  the  HE  case  the  fraction  is  initially  about 
one-half  and  decreases  monotonically  thereafter  at  about  the 
same  rate  as  in  the  nuclear  case. 


Section  5 
CONCLUSIONS 


We  have  described  numerical  simulations  carried  out 
for  a  600-ton  HE  burst  and  a  1-kton  nuclear  burst,  both  at 
166  ft.  The  code,  gridding,  and  method  of  solution  are  the 
same  in  the  two  calculations.  Although  the  shock  waves  in 
the  two  cases  propagate  at  roughly  the  same  speed  and  break 
away  at  roughly  the  same  time,  and  although  the  pressure 
fields  relax  to  ambient  in  similar  fashion,  we  find  signifi¬ 
cant  differences,  of  which  the  following  appear  to  be  the 
most  important. 

(i)  The  HE  flow  velocities  are  systematically 
larger. 

(ii)  In  the  regular  reflection  region  (underneath 
and  close  to  the  fireball),  the  HE  case 
exhibits  two  overpressure  peaks  at  the 
surface,  rather  than  one,  due  to  re-reflec¬ 
tion  of  the  reflected  shock  from  the  contact 
surface  between  air  and  detonation  products.- 

(iii)  For  the  HE  case  the  upper  vortex  forms  first, 
followed  by  the  reverse  vortex  near  the  axis 
of  symmetry  and  the  ground.  Adjacent  HE 
products  begin  to  be  entrained  into  a 
positive  vertex  over  a  longer  period  of  time, 
several  seconds.  In  the  nuclear  case  the 
negative  vortex  is  the  dominant  one;  it  is 
larger  then  the  HE,  persists  longer,  and 
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contains  larger  velocities  than  the  positive 
vortex  at  comparable  times. 

(iv)  The  HE  flow  establishes  a  pattern  of  four 
vortices,  two  forward  and  two  reversed, 
instead  of  one  of  each. 

(v)  The  stem  of  the  nuclear  fireball  appears  to 
exhibit  a  Rayleigh-Tayior  instability,  absent 
in  the  HE  case. 

Since  the  velocities  on  axis  are  higher  in  the  HE 
case  (the  upper  positive  vortex  is  larger),  fireball  rise  is 
faster  than  in  the  nuclear  case.  The  nuclear  case,  however, 
probably  scours  up  more  dust  because  the  velocities  are 
larger,  and  the  reverse  vortex  is  larger  and  more  persistent. 
It  is  difficult  to  argue  conclusively  on  this  point  because 
so  much  depends  on  terrain,  conditions  in  the  boundary  layer, 
and  other  physical  effects  not  included  in  this  model  (e.g., 
precursor  heating  and  turbulence).  Further  study  of  the 
tracer  particle  motions  we  have  calculated  is  expected  to  be 
illuminating  in  this  regard. 
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ADAPTATION  CF  FLUX-CORRECTED  TRANSPORT  ALGORITHMS  FOR  WCDELING  XSTT  FLOWS 

M.  A.  Fry,  Science  Applications,  Inc. 
and 

D.  L.  Book,  Naval  Research  Laboratory 

INTRODUCTION 

In  this  paper  we  describe  a  series  of  calculations  carried  out  as  part  of  an 
ongoing  effort  aimed  at  studying  blast  wave  diffraction  effects  in  air.  The 
phenomena  of  chief  interest  to  us  include  velocity  fields,  particularly  those 
associated  with  the  toruses  (both  forward  and  reverse)  in  the  neighborhood  of  the 
rising  fireball,  and  the  distribution  of  dust  lifted  off  the  ground  by  the  winds 
and  the  structure  of  the  cloud  at  the  time  of  stabilization,  we  are  interested  in 
studying  the  nature  of  the  gas-dynamic  discontinuities  which  appear,  the  vortices 
(both  forward  and  reverse),  and  how  the  dust  content  of  the  air  affects  the 
evolution  of  the  blast  wave. 

The  technique  we  have  employed  for  this  purpose  is  numerical  mode  line .  One- 
and  two-fluid  hydrocodes  based  on  the  Flux-Corrected  Transport  (PCT)1  shock¬ 
capturing  techniques  have  been  used  to  simulate  airblast  phenomena  in  one  and  two 
dimensions.  FCT  refers  to  a  class  of  state-of-the-art  fluid  computational 
algorithms  developed  at  NRL  in  the  course  of  the  past  ten  years  with  supersonic 
gas-dynamic  applications  expressly  in  mind.  We  have  concentrated  on  model ina  the 
"Direct  Course"  event,  an  experiment  to  be  fielded  shortly  by  the  Defense  Agency: 
a  600-ton  ammonium  nitrate  +  fuel  oil  (ANFO)  charge  is  detonated  at  a  height  of 
burst  (BOB)  of  166  ft.  The  results  are  most  conveniently  exhibited  in  terms  of 
velocity  vector  plots  and  tracer  particle  trajectories.  Examples  of  these  are 
presented  to  illustrate  our  results  and  conclusions. 

The  plan  of  the  paper  is  as  follows:  In  the  next  section  we  discuss  our 
numerical  techniques  and  validation  procedures.  In  Section  3  we  describe  the  600- 
ton  2D  KB  calculations.  In  Section  4  we  summarize  our  conclusions  and  discuss 
their  domain  of  validity. 

Manuscript  approved  September  13,  1983. 
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NUMERICAL  TREATMENT 


As  described  by  Boris  and  Book1 ,  PUT  is  a  finite-difference  technique  for 
solving  the  fluid  equations  in  problems  vhere  sharp  discontinuities  arise  (e.a. , 
shocks,  slip  surfaces  and  contact  surfaces).  It  modifies  the  linear  properties  of 
a  second-  (or  higher-)  order  algorithm  by  adding  a  diffusion  term  during  convec¬ 
tive  transport,  and  then  subtracting  it  out  "almost  everywhere"  in  the  antidiffu¬ 
sion  phase  of  each  time  step.  The  residual  diffusion  is  just  large  enough  to 
prevent  dispersive  ripples  from  arising  at  the  discontinuity,  thus  ensuring  that 
all  physically  positive  conserved  quantities  retain  positive.  FCT  captures  shocks 
accurately  over  a  wide  range  of  parameters.  No  information  about  the  number  or 
nature  of  the  surfaces  of  discontinuity  need  be  provided  prior  to  initiatina  the 
•calculation. 

The  FCT  routine  used  in  the  present  calculations,  called  JPBFCT  (an 
advanced  version  of  ETBPCT2),  consists  of  a  flexible,  general  transport  module 
which  solves  1-D  fluid  equations  in  Cartesian,  cylindrical,  or  spherical  geometry. 
It  provides  a  finite-difference  approximation  to  conservation  laws  in  the  general 
form: 


•iL.  /  ?dv  *  -  /  $>(u-u  )*dA  +  /  tdA 

\  £  n  /  2  *  t  a.  \ 


«V(t) 


«A(t) 


«A(t) 


where  $  represents  the  mass,  momentum,  energy  or  species  mass  density  in  cell 
oV(t),  ju  and  Ug  represent  the  fluid  and  grid  velocities,  respectively,  and  t 
represents  the  pressure/vork  terms.  This  formulation  allows  the  grid  to  slide 
with  respect  to  the  fluid  without  introducing  any  additional  numerical  diffusion. 
Thus,  knowing  where  the  features  of  greatest  interest  are  located,  one  can  con¬ 
centrate  fine  zones  where  they  will  resolve  these  features  most  effectively  as  the 
system  evolves. 


The  same  transport  routine  was  employed  for  both  coordinate  directions  In  the 
2D  r-z  code  (called  FAST2D)  via  timestep  splitting.  A  Jones-Wilkins-Lee  (JWL) 
equation  of  state  (EOS)  was  used  for  the  detonation  products  and  a  real -air  EDS 
was  used  outside  the  HE-air  interface3 .  The  routine  was  written  in  the  form  of  a 
cable  lookup,  using  interpolation  with  logarithms  to  the  base  16  computed  by  means 
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of  logical  snifts1* .  By  thus  taking  account  of  the  architecture  of  the  machine  :  in 
these  calculations,  a  32-bit -word  two-?ipe  Ttexas  Instruments  AS C)  it  was  possiole 
to  generate  very  efficient  vector  code,  decreasing  the  time  required  for  EDS  cal¬ 
culations  to  a  snail  fraction  of  that  required  for  the  hydro.  The  EDS  specifies 
pressure  as  a  function  of  density  and  internal  energy.  In  mixed  cells  the 
combined  pressure  was  calculated  according  to  Dalton's  law. 

The  initial  conditions  were  taken  to  be  the  self-similar  flow  field  used  by 
Kuhl,  et  al.~ ,  corresponding  to  a  spherical  Chapman-Jououet  detonation  at  the  time 
the  detonation  wave  reaches  the  charge  radius  (Fig.  1).  This  was  prooaaated  with 
the  ID  spherical  code  until  the  detonation  front  attained  a  radius  just  smaller 
than  the  HOB,  at  which  tine  the  solution  was  laid  dc*n  on  the  2D  mesh  (Fig.  2). 

The  boundary  conditions  were  cdiosen  to  enforce  perfect  reflection  on  the 
ground  and  on  the  axis  of  symmetry  £(<3*/^)^  *  where  <t>  *  P,o,vt,  and  v^  <=  0] , 
where  "t"  and  "n"  denote  tangential  and  normal  components,  respectively,  and 
outflow  on  the  outer  and  the  top  boundaries  [(d$/dn)^  »  0,  where  <f>  *  o ,p,vt,vn\ 

For  the  2D  calculations  the  mesh  was  typically  200  x  100.  Fixed  qriddina  was 
used  to  minimize  numerical  errors.  The  zone  sizes  were  2.1  m  x  2.1  m,  respec¬ 
tively.  Ftor  the  late-time  calculations,  a  fixed  mesh  with  100  zones  in  the  radial 
and  200  zones  in  the  vertical  direction  was  used,  with  all  cells  of  dimension 
4.2  m  x  4.2  m. 

Tto  study  the  notion  of  dust  particles  in  the  flow  field  generated  by  the 
calculation,  the  simplest  model  describes  dusty  air  as  a  single  phase  with  density 
and  adiabatic  index  chosen  appropriately.  This  approach  ignores  the  properties 
associated  with  the  particulate  structure  of  the  dust  and  the  process  of  scouring 
by  which  the  dust  enters  the  air.  A  more  realistic  picture  results  if  we  treat 
the  dust  as  a  distinct  phase,  described  by  equations  of  mass,  momentum  and  energy 
conservation,  as  has  been  done  in  one  dimension  by  Miura  and  Glass5.  The  dust 
equations  are  aoupled  to  those  describing  the  air  through  drag  and  heat  transfer 
terms. 


) 


It  has  been  pointed  out  by  Kuhl  that  in  such  a  treatment  a  dust  oar- 
ticle  tends  to  become  entrained  in  the  prevailing  flow  over  a  distance  ~  10 5 
particle  diameters.  Thus  a  one-ohase  description  is  satisfactory  whenever  oar- 
ticle  sizes  are  at  least  a  thousand  times  snaller  than  the  snallest  lenoth  scale 
in  the  hydrodynamics.  For  the  present  calculation,  this  scale  is  roughly  1  m,  so 
particles  snaller  than  1  rim  can  be  regarded  as  totally  entrained. 

When  the  mass  density  of  the  dust  component  is  snail  comoared  with  air 
(or  HE  product)  density,  a  further  simplification  is  possible.  Entrained  dust 
particles  can  be  followed  by  passive  advection.  That  is,  the  wind  fields  ux, 
Uy  are  taken  from  a  dust-free  hydrodynamics  calculation,  and  dust  is  advected  in 
these  fields  according  to 


In  this  approximation  we  ignore  the  effect  of  the  momentum  and  energy  transfer  on 
the  air  phase.  The  same  approximation  can  be  used  for  larger  (nonentrained) 
particles  also,  provided  we  include  inertia  and  drag  by  using  the  force  law 


3t  *  "  "SSz  +  D(-  " 


where  V  is  the  particle  velocity,  g  is  the  acceleration  due  to  gravity,  and  D  is 
the  enpirical  drag  coefficient  employed  by  Miura  and  Glass5. 

Equations  (2)  and  (3)  apply  best  in  the  limits  of  extremely  snail  and 
extremely  large  particles,  respectively.  Although  they  restrict  the  scope  of  the 
calculation  (by  requiring  the  dust  content  to  be  small ) ,  they  have  the  computa¬ 
tional  advantage  of  allowing  us  to  obtain  time-dependent  dust  distributions  for 
many  different  choices  of  dust  size  spectrum  and  initial  distribution  from  a 
single  hydrodynamics  calculation. 
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600-7TN  ANTO  EXPLOSION  AT  166  FT 


The  yield  and  HOB  (600  tons  and  166  ft,  respectively)  in  the  calculation  were 
chosen  to  equal  the  values  used  in  the  Defense  Nuclear  Aaency  Direct  Course 
experiment,  which  ve  are  simulating.  The  Chapman -Jouguet  parameters  used  to 
initialize  the  spherical  free-field  calculation  were  taken  to  be  those  for  the 
NHuND,-fuel  oil  (ANPO)  mixture  used  as  the  explosive. 

Figs.  3(a)-(c)  show  the  contours  of  HE  density  and  internal  energy  per  unit 
mass  and  the  velocity  arrow  plot  at  t*0,  just  before  the  reflection  at  ground  zero 
occurs.  Figs.  3(d)-(f)  show  the  corresponding  plots  54  ms  later,  while  Pigs. 
3(g)— (i)  show  them  after  245  ms.  Note  the  reflected  shock  proceeding  uoward, 
reflecting  again  off  the  fireball,  and  propagating  back  in  a  downward  and  outward 
direction.  The  interaction  of  this  shock  with  the  radially  inward  flow  near  the 
ground  generates  the  reverse  vortex,  vhich  is  clearly  seen  in  Fig.  3(i).  Note 
also  the  positive  vortex  forming  near  the  top  of  the  grid  in  the  same  olot.  The 
latter  results  when  the  upward-propagating  reflected  shock  interacts  with  the 
radially  outward  flow  near  the  top  of  the  fireball?  it  is  not  produced  by  the 
buoyant  rise  of  the  fireball,  vrtuch  at  these  early  times  has  scarcely  begun. 

To  look  at  the  evolution  of  the  fireball  at  late  times,  we  reinitialized  on  a 
larger,  coarser  grid,  representing  a  cylinder  400  m  in  radius  and  800  m  high.  The 
first  300  cycles  approximately  reproduce  the  early-time  results.  The  spherical 
shock  breaks  away  and  leaves  the  mesh.  The  flows  remaining  on  the  grid  are  now 
subsonic  everywhere.  Then  the  fireball  begins  to  rise  and  the  subsequent  develop¬ 
ment  is  due  to  the  combination  of  buoyant  rise  and  the  action  of  the  vortices  set 
up  by  the  early  shocks. 

Figure  4a,  which  shews  the  trajectories  of  passively  advected  tracer  oar- 
ticles  <wer  the  time  interval  1.8  sec  to  3.97  sec,  displays  the  vortices  very 
clearly.  Fig.  4b  shows  the  particle  paths  for  the  time  interval  3.97  sec  to  7.34 
sec.  Notice  that  there  are  four  vortices  visible  in  the  plot:  two  positive  and 
two  reversed.  The  additional  small  vortices  are  apparently  a  consequence  of 
entrainment  by  the  major  ones.  As  far  as  we  know,  their  existence  has  not  been 
noted  previously. 
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It  is  clear  that  the  major  qualitative  features  persist  longest  in  the 
velocity  plots.  This  is  not  surprising,  as  the  circulation  patterns  represented 
by  tne  vortices  have  essentially  infinite  lifetimes  in  the  absence  of  viscosity. 
We  have  run  out  to  stabilization  (not  shown  here)  and  have  found  that  these 
features  persist  in  the  velocity  plots  to  the  very  end.  At  all  times  t>0  the  oeak 
flow  velocity  in  tne  HE  case  exceeds  that  in  the  comparable  point  source  solution. 
This  is  a  reflection  of  the  fact  that  the  Chapman -Joug ue t  solution  at  a  radius  of 
10  m  has  a  pressure  peak  of  52  tear,  vs  3  tear  for  the  Sedov  solution  at  the  same 
radius.  The  means  that  the  former  starts  out  with  much  more  violent  motion,  i.e. , 
fluid  velocities  an  order  of  magnitude  larger.  All  in  all,  in  many  respects  the 
HE  case  does  not  closely  resemble  a  point  source. 


CONCLUSIONS 

We  have  described  a  numerical  simulation  of  the  Direct  Course  Event.  The 
code,  gridding,  and  method  of  solution  are  the  same  in  the  two  calculations.  The 
following  conclusions  appear  to  be  among  the  most  important. 

(i)  The  flow  establishes  a  pattern  of  four  vortices,  two  forward  and  two 
reversed,  instead  of  one  of  each. 

(ii)  The  upper  vortex  forms  first,  followed  by  the  reverse  vortex  near  the 
axis  of  sywnetry  and  the  ground.  Adjacent  HE  products  beain  to  be 
entrained  into  a  positive  vortex  over  a  longer  period  of  time,  several 
seconds. 

(iii)  Cnee  picked  up  (scoured)  off  the  ground,  dust  is  efficiently  trans¬ 
ported  upward  by  the  reverse  vortex  farther  fran  the  axis. 

Phenomena  neglected  in  the  present  model  (e.g.,  terrain,  conditions  in  the 
boundary  layer,  turbulence,  humidity,  etc.)  are  unlikely  to  alter  the  above  con¬ 
clusions,  vdiich  mainly  depend  cn  the  characteristics  of  the  solutions  in  the 
interior  of  the  mesh  and  over  long  periods  of  time.  Further  study  of  the  tracer 
particle  motions  we  have  calculated  is  likely  to  be  illuminating,  particularly 
when  we  begin  to  consider  the  evolution  of  various  initial  configurations  as  a 
function  of  particle  size.  In  closing,  it  is  appropriate  to  emphasize  the  far- 
reaching  significance  of  the  role  played  by  the  HE-air  interface  in  the  dynamics 
of  both  airblast  and  cloud  rise  phenomena,  and  the  importance  for  numerical 
simulation  of  correctly  treating  this  interface. 
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The  self-similar  solutions  to  the  problem  of  a  planar  shock  with  Mach  number  M$  reflecting 
obliquely  from  a  wedge  with  vertex  angle  are  obtained  to  arbitrary  accuracy  by  expanding  the 
Quid  quantities  in  power  series  in  the  scaled  variables  {  *  x/t,  tj  -  y/t.  For  single  Mach  reflection, 
there  are  four  distinct  regions:  (a)  the  ambient  gas  ahead  of  the  incident  shock;  (b)  the  gas  behind 
the  incident  shock  and  outside  the  reflected  (bow)  shock;  (c)  the  region  bounded  by  the  Mach  stem, 
the  wedge,  and  the  contact  surface  (slip  line)  extending  from  the  triple  point;  and  (d)  the  doubly 
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20.  ABSTRACT  (Continued) 

shocked  medium  bounded  by  the  contact  surface,  wedge,  and  reflected  shock.  In  region 
(b)  the  solution  is  known  immediately  in  terms  of  Mv  6  w,  and  the  conditions  in  (a). 

Hie  shapes  of  the  Mach  stem,  contact  surface  and  bow  shock  are  expressed  parametrically 
as  |  -  F(s),  r?  -  G(s).  Then  p^  uc,  pe,  and  pd,  ud,  vd,  pd  are  obtained  by  expanding 
variables  in  double  power  series.  e.g., 

pcd,  i ?)  - 1 

ij  i] 

substituting  in  the  ideal  fluid  equations,  and  equating  coefficients  of  like  powers  through 
some  order  N  ■  max(H).  The  resulting  algebraic  equations  are  solved  subject  to  the 
additional  relations  obtained  by  applying  the  reflection  conditions  on  the  wedge,  together 
with  the  jump  conditions  on  the  boundaries  ac  and  bd,  approximated  by  power  series 
expansions  of  die  F  and  G  functions.  Since  all  these  equations  are  nonlinear,  solutions 
are  obtained  by  iteration  with  N  increasing  until  convergence  is  obtained.  The  Ben-Dor 
equation  for  the  fluid  quantities  in  regions  c,  d  at  the  triple  point  is  used  to  give  initial 
values.  Because  variation  within  each  region  is  smooth,  effectively  exact  descriptions  of 
most  features  of  interest  can  be  obtained  using  series  with  £  20  terms.  There  are  thus 
<  200  quantities  in  the  discretization  of  the  problem,  compared  with  ~  104  in  a  con¬ 
ventional  finite-difference  treatment.  The  method  generalizes  readily  to  complex  and 
double  Mach  reflections. 
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INTRODUCTION 

ZHie  problem  of  a  planar  shock  reflecting  from  an  oblique  surface  goes 
back  over  a  hundred  years  to  Ernst  Mach.  Although  this  problem  is  important 
in  its  own  right,  much  of  the  interest  in  it  arises  because  of  the  need  for 
better  understanding  of  Mach  reflection  in  more  complicated  situations.  The 
field  has  been  the  object  of  particular  interest  during  the  last  thirty 
years;  the  experimental  and  theoretical  research  carried  out  during  this 
period  have  been  reviewed  by  Ben-Dor.1 

A  constant  planar  shock  propagating  into  a  uniform  ambient  gas  gives 
rise  in  the  absence  of  reflection  to  a  second  medium  with  uniform  thermody¬ 
namic  properties  in  the  region  behind  the  shock.  If  it  propagates  in  a  shock 
Manuscript  approved  September  13,  1983. 


tube  whose  walls  are  not  parallel  to  the  di recoion  of  propagation  (because  a 
wedge  has  been  inserted  along  the  side) ,  a  reflected  shock  wave  propacates 
back  into  the  interior  of  the  shock  tube.  When  the  wedge  angle  is  large,  so 
that  the  primary  shock  is  incident  nearly  normally  on  it,  the  reflected  shock 
is  also  planar  and  no  other  gasdynamic  discontinuities  appear  near  the 
reflection  point.  As  the  'wedge  angle  decreases,  so  that  the  incident  3hock 
becomes  more  and  more  nearly  glancing,  it  becomes  impossible  for  a  fluid 
particle  to  traverse  both  incident  and  reflected  shocks  and  still  "turn  the 
corner"  enough  to  end  up  moving  parallel  to  the  wedge  surface.  At  least  one 
additional  shock  (the  Mach  stem)  must  appear  (Fig.  1),  intersecting  the 
others  at  a  so-called  triple  point.  3ecause  some  of  the  material  in-  the  zone 
between  the  Mach  stem  and  the  reflected  wave  has  been  shocked  once  and  some 
twice,  another  gasdynamic  discontinuity  (a  contact  surface)  must  also  extend 
from  the  triple  point  into  this  region,  terminating  somewhere  on  the  wedge 
surface,  the  reflected  shock  may  terminate  at  the  corner  of  the  wedge 
(attached  shock),  or  upstream  from  this  point  (detached  shock),  or  may  run 
into  a  second  triple  point  between  the  first  one  and  the  corner  (double  Mach 
reflection) .  The  latter  case  occurs  in  general  for  smaller  wedge  angles  than 
does  single  Mach  reflection;  an  intermediate  case  (complex  Mach  reflection) 
is  also  observed. 

Cne  would  like  to  derive  a  theoretical  description  of  Mach  reflection 
which  would  complement  the  experimental  results  and  address  some  of  the 
questions  the  latter  leave  unanswered,  such  as  the  structure  of  the  contact 
surface  near  the  wedge  surface,  and  whether  a  "triple  Mach"  regime  exists. 

An  analytical  solution  is  out  of  the  question,  though  pieces  of  the  problem 
(e.g.,  the  flow  in  the  neighborhood  of  the  triple  point*)  can  be  solved. 
Recourse  must  therefore  be  had  to  numerical  methods. 
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Straightforward  differencing  of  the  fluid  equations  (using  either  an 
ideal  or  realistic  equation  of  state)  has  been  carried  out  with  considerable 
success*'  .  A  rectilinear  Sulerian  mesh  with  ~i  00  -  200  zones  along  each 
axis,  possibly  varying  in  size  as  a  function  of  position  to  improve  resolu¬ 
tion  near  the  Mach  stem,  is  used.  The  conditions  ahead  of  and  behind  the 
incident  shock  are  used  as  boundary  conditions,  along  with  a  reflection  con¬ 
dition  on  the  wedge.  State-of-the  art  shock  capturing  techniques  such  as 
Flux-Corrected  Transport  (FCT)3  resolve  the  envelope  (reflected)  waveform 
accurately  and  permit  almost  all  the  other  gasdynamic  diseontinuties  in  the 
problem  to  be  distinguished  (Fig.  2).  To  date,  however,  such  numerical 
solutions  have  not  surpassed  experimental  interferometric  data  in  accuracy, 
nor  have  they  succeeded  in  answering  any. of  the  outstanding  questions 
associated  with  Mach  reflection. 

Moreover,  such  calculations  leave  a  distinct  impression  of  brute  force. 
Advancing  the  four  fluid  equation  103  -  104  timesteps  on  10“  or  more  zones 
seems  profligate,  particularly,  in  view  of  the  smoothness  of  most  of  the  gas- 
dynamic  discontinuities  and  the  gentle  variation  of  the  fluid  quantities  in 
the  regions  they  bound.  In  much  of  the  domain  of  the  calculation  the  solu¬ 
tion  is  known  a_  priori  and  does  not  change  in  time . 

Furthermore,  the  desired  solution  is  actually  self -similar .  In  the 
scaled  variables  Z  *  x/t,  n  »  y/t,  where  t  is  time  and  the  origin  of  the 
coordinate  system  is  fixed  at  the  wedge  comer,  the  observed  wave  forms  are 
stationary.  It  is  thus  natural  to  seek  a  solution  to  the  ideal  fluid  equa¬ 
tions  in  terms  of  these  similarity  variables. 


In  the  present  paper  we  treat  the  shock -on -wedge  problem  by  expanding 
the  fluid  variables  and  the  functional  forms  of  the  boundaries  in  power 
series  in  £  and  n.  The  expansion  coefficients  are  found  by  imposing  the 
boundary  conditions  (Rankine-Hugoniot  conditions  on  the  shock,  perfect 
reflection  on  the  wedge  surface )  . 

Mumerous  fluid  dynamics  problems  have  been  solved  by  power  series  tech¬ 
niques,  as  described  by  Van  Dyke.4  Here  we  follow  the  general  approach  he 
advocates:  first  do  the  lowest-order  term s  "by  hand,"  then  afterward  auto¬ 

mate  the  procedure,  finally  using  sophisticated  mathematical  techniques  to 
find  the  limiting  values  of  critical  numbers  (e.g.,  boundaries  between  two 
reflection  regimes).  Unlike  most  of  the  problems  described  by  Van  Dyke,  the 
algebraic  equations  determining  the  power  series  coefficients  here  are  highly 
nonlinear  and  need  to  be  solved  by  iteration;  the  proper  iteration  scheme  is 
not  obvious,  but  must  be  found  by  experimentation. 

In  the  following  section  we  derive  the  equations  to  be  solved.  The  next 
section  describes  the  iteration  procedure,  the  program  implementing  it,  and 
the  display.  The  final  section  summarizes  our  conclusions  and  discusses 
extension  of  the  calculation  to  higher  order  and  to  the  more  complicated  Mach 
reflection  cases . 

DERIVATION  OF  EQUATIONS 

In  what  follows  we  use  the  label  "a"  for  quantities  in  the  ambient  gas 
(ahead  of  the  incident  shock);  "b"  behind  the  incident  shock;  "c"  between  the 
contact  surface  and  the  Mach  stem;  and  "d"  for  the  doubly  shocked  region 
between  the  reflected  shock  and  the  contact  surface.  Surfaces  of  gasdynamic 
discontinuities  are  labelled  by  the  two  regions  they  separate,  e.g.  "cd"  for 
the  contact  surface.  The  triple  point  is  labeled  by  the  subscript  "T, "  and 
the  origin  is  at  the  point  of  attachment  (the  wedge  corner) . 


In  terms  of  the  similar! ty  variables  i,  n,  the  fluid  equations  become 
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where  y  is  the  adiabatic  index. 


we  define  dimensionless  variables  bv 


X  »  ;/c,  Y  *  n/c  ; 


d  *  (u-S)c*1  ,  v  «  (v-n)c”1  , 


p/p a  f  p  -  P/Pa  f 


where 


c  “  Vpa' 


then  the  fluid  eouations  become 
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In  each  of  regions  c  and  d  we  expand  R,  0,  V,  P  in  double  power  series, 


R(X,Y)  -  P00  +  R10X  +  R0iY  +  R20X2  *  R,,XY  +  Rn2Y2  - 


*  r30x3  +  r2,x2y  *  r,2xy2  ♦  r03y3  +  .  .  .  . 


U(X,Y)  -  a00  +  u10x  +  O01Y  *  d20x2  +  uuxy  *  o02y2  + 


+  0anx3  *  n,,X2Y  *  U,,XY2  ♦  0n,Y3  +  .  .  .  . 


V( X , '£)  -  V00  +  V10X  +  V01Y  +  V20XZ  ♦  VnXY  +  VQ2Y2  + 


'T  v  3 


♦  V,„XJ  +  V,,X‘Y  +  V,  „XY‘  +  Vn,YJ  +  .  .  .  . 
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P(X,Y)  -  P00  +  ?,0X  +  P31Y  +  P20XZ  +  PUXY  +  P02Y2  + 


p30x3  *  ?21x2y  ♦  ?.2xy;  ♦  ?03x3  ♦  .  .  .  . 


To  apply  the  reflection  boundary  condition,  we  set 
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It  follows  that  we  must  have 


R01  *  R11  *  R 21  *  R  31 


•  •  .  *  0  ; 


^01  *  Uli  =  u2i  *  "si  *  *  *  *  *  O' 


?01  *  P11  *  P21  "  ?31  *  *  *  *  =  0; 


V00  *  V10  *  V20  *  v30 . °‘ 


To  proceed  further,  we  substitute  Eqs.  (13)-{16)  in  Eqs.  (9)— (12)  and 
collect  terms  of  like  powers  X1  ,  equating  their  coefficients  to  zero  order 
by  order.  The  task  of  doing  this  by  conventional  methods  rapidly  becomes 
hopeless,  but  it  can  readily  be  automated.  To  see  what  is  going  on  in  the 
lowest  orders  of  our  "hand"  calculation,  we  expand  using  MACSYMA,  the  sym¬ 
bolic  manipulation  program  developed  at  the  Math lab  of  MIT.  This  works  up  to 
about  order  N  «  6,  where  we  define  the  order  of  the  expansion  by  N  « 
max(i+j ) .  Beyond  that  point  storage  limitations  make  it  necessary  to  use 
"tricks,"  and  eventually  terminate  the  process  entirely. 


When  this  is  done,  certain  patterns  emerge.  We  find  that  we  must  have 
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Some  of  the  equations  then  vanish  identically.  In  each  of  regions  c  and  d 


there  are  more  variables  than  nontrivial  equations,  as  shown  by  the  following 


table: 


Eouations 


variables 


Order 


0  V 


Total 


Total 


1  1  1 


Table  1 


Table  1  shows  that  in  each  region  there  are  three  more  unknowns  than 


equations  in  each  order,  except  for  N  «  0,  where  the  number  is  four.  The 


additional  equations  are  supplied  by  the  jump  conditions,  imposed  on  surfaces 


ac  and  bd  for  regions  c  and  d,  respectively. 


If  we  represent  a  shock  boundary  parametrically  by  x»  F(s),  Y  •  S(s), 


then  the  Rankine-Huooniot  conditions  become 


X*  (U  -  U)  +  Y'  (V  -  V)  »  0; 


Y'  (R  0  -  R  C)  *  X'  (R  V  -  R  V); 


?  +  R  O2  +  V2)  «  ?  -  R  CJ2  -  V2); 


'  26) 


U2 


V" 


21?  =  ~2  *  -2  +  2YP 

Y-1  y_T 


while  on  a  contact  surface 


P  -  (28) 

Y'  U  -  X'  V;  (29) 

Y'  0  «  X’  V.  (30) 

Here  the  variables  in. front  of  and  behind  the  discontinuity  are  denoted 
by  unbarred  and  barred  symbols,  respectively. 

we  know  that  ac  is  nearly  a  straight  line  normal  to  the  wall,  so  the 
representation 


X  -  XQ  +  XlY  +  X2Y2  +  X3Y3  +  .  .  .  .  (31) 

a  a 

is  valid.  Substituting  Eg.  (31)  and  the  ambient  conditions  u  *  v  *  0  in 
Eg.  (24),  the  condition  for  continuity  of  parallel  velocities  at  ac,  which 


takes  the  form 


[I] 


(V01  +  15  Y1  +  U1 C  +  ‘ 


-  1)YiY2  *  (Vn  +  Ua02) 


Vj  3  &  b  b 

etc.  (Note  that  u  and  v  ,  like  u  and  v  ,  are  constants,  whereas  tJ  ,  V  , 


ua,  Va  are  not.)  Substitution  of  35.  (35)  is  Bqs.  (25)— (27)  with  X'  »  1  then 


yields  three  conditions  per  order  on  the  "d"  variables. 


At  this  point  we  need  three  more  conditions:  two  for  the  remaining 


variables  in  order  zero,  and  another  to  determine  Xn.  For  this  purpose  we 


use  Sqs.  ( 23) -( 30) .  Note  that  to  apply  these  conditions  everywhere  on  cd 


would  overdetermine  the  system.  This  is  equivalent  to  asserting  the 


impossibility  of  continuing  the  solution  across  ab  from  a  to  b,  across  ac 


from  a  to  c,  across  bd  from  b  to  d,  and  then  across  cd  from  c  back  to  d. 


where  it  is  already  known.  Yhe  only  freedom  we  have  is  to  impose  Sqs.  (28)- 


(30)  at  the  triple  point. 


A  il.** V./. 


r 


deteraines  the  ratio 


For  our  trial  calculation  we  work  to  order  N  ■  5.  Inclucina  the  cuanti 


ties  which  vanish  identically,  we  then  have  163  equations  in  183  unknowns, 


which  are  tc  be  solved  simultaneously. 


SOLUTION  OF  EQUATIONS 


The  solution  is  obtained  by  iteration.  Several  points  are  important  in 


the  design  of  a  satisfactory  iteration  scheme: 


(i)  Good  values  of  the  quantities  X0,  X*,,  YT,  etc.,  can  be 


obtained  from  the  experimental  data  (Fig.  1)  and  used  as  initial 


eraesses , 


(ii)  The  system  can  be  made  quasilinear  if  we  solve  order  by  order. 


starting  with  N  ■  0,  and  using  in  any  given  order  the  previously  obtained 


values  of  all  variables  not  being  solved  for  in  that  order. 


(iii)  To  reduce  the  effect  of  possible  instability  in  parts  of  the 


scheme,  all  variables  are  updated  using  some  form  of  (under Relaxation,  e.c. 


new  value  ■  old  value  +  relaxation  factor  *  corrections. 


The  current  version  of  the  code  works  as  follows: 


(i)  For  N  m  0,  use  the  X-  and  Y-  independent  terms  of  Eqs.  (9) -(12), 


together  with  the  zeroth-order  form  of  Eqs.  ( 25 ) — ( 27 )  and  Bq.  (43)  to  update 


c  c  c  c  c  c  c 
R  ,  R  ,  U  ,  U  ,  V  ,  P  ,  P  ; 

00  10  00  10  01  00  10 


(ii)  Do  the  same  thing  [using  Bq.  (44)  instead  of  (43)]  for  R  , 


.,  Pd  ; 
10 


(iii)  Solve  Bqs.  (9) — ( 1 2)  and  Bqs.  (25)-(27)  on  ac  order  by  order  for 


N  »  1  to  5  to  obtain  the  region-c  variables; 


(iv)  Solve  Bqs.  (9)— (12)  and  Bqs.  (2*'-(27)  on  bd  order  by  order  for 


N  ■  1  to  5  to  obtain  the  region-d  variables  plus  Y^.,  ; 


(v)  Iterate  Bqs.  (40) -(42)  together  with  the  successive  orders  of 


Bq.  (32)  [e.g.,  Bq.  (34)  for  X2]  to  solve  for  XT,  Y^,  and  X,. ,  i  ■  0,  1 


6.  (To  a  good  approximation,  X^  »  0  for  all  i  *  0,  2.) 


This  scheme  is  repeated  until  no  further  change  (to  some  preset  tolerance)  in 
the  variables  occurs.  The  program,  written  in  Fortran,  runs  on  a  VAX  11/780 
at  about  one  second  per  iteration. 

Results  are  most  conveniently  displayed  as  contour  plots  in  pressure  (or 
density).  Although  it  is  possible  to  drive  the  plotting  pen  directly  using 
the  exact  formula  ?(X,Y),  it  is  easier  to  declare  a  very  large  array  (e.g., 
500*500),  fill  it  with  pressures  calculated  at  every  X,  Y,  and  use  a  standard 
contour  plotting  package  on  this. 

CONCLUSION 

The  program  described  above  is  still  under  development,  and  no  results 
have  yet  been  generated.  For  this  sample  problem,  it  seems  straightforward 
to  carry  the  method  to  a  successful  solution.  We  have  encountered,  and 
apparently  overcome,  two  types  of  difficulties:  determining  the  formulation 
for  the  problem,  and  solving  the  resulting  set  of  equations.  Only  if  both 
stages  are  handled  correctly  will  useful  answers  results.  There  remains  the 
(programming)  task  of  automating  the  solution,  so  as  to  work  to  arbitrarily 
high  order  N.  This  is  of  interest  chiefly  in  connection  with  studying  the 
behavior  of  the  roll-up  in  the  contact  surface. 

Finally,  extension  to  other  forms  of  Mach  reflection  is  of  interest. 
Attached  double  (or  triple)  Mach  reflection  presents  no  problem  in  principle, 
and  can  be  handled  using  the  techniques  described  here.  Complex  Mach 
reflection  and  detached  shocks  are  less  clear.  At  present  we  do  not  know  how 
to  formulate  the  problem  in  either  situation-  so  as  to  make  it  well -posed. 

This  will  be  addressed  in  future  -work. 
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